Set Up

library(tidyverse)
library(terra)
library(sf)
library(mapview)
library(raster)
library(rgeos)
library(lubridate)
library(ggplot2)
library(exactextractr)
library(patchwoRk)
library(gridExtra)

Define Fire Parameters

USDA National Forest Type Group Dataset (2004)

Conifer Forest Type Groups: Douglas-Fir, Ponderosa Pine, Fir-Spruce-Mountain Hemlock, Lodgepole Pine

# forest type groups and key
conus_forestgroup <- raster('data/conus_forestgroup.tif')
forest_codes <- read_csv('data/forestgroupcodes.csv')

# set crs
crs = crs(conus_forestgroup)

EPA level-3 Ecoregions

Canadian Rockies, Idaho Batholith, Middle Rockies, Columbian Mountains - Northern Rockies

# level 3 ecoregions
l3eco <- st_read('data/us_eco_l3.shp') %>% 
  st_transform(., crs=crs)
## Reading layer `us_eco_l3' from data source 
##   `G:\Other computers\My Laptop\Documents\Grad School\Research\ConiferRegeneration\data\us_eco_l3.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1250 features and 13 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2356069 ymin: 272048.5 xmax: 2258225 ymax: 3172577
## Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic_USGS_version
# select northern rocky mountains from level3 ecoregions
eco_select <- l3eco %>% 
  filter(NA_L3NAME %in% c('Canadian Rockies','Columbia Mountains/Northern Rockies','Middle Rockies','Idaho Batholith'))

mapview(eco_select)

MTBS Fires Boundaries

Criteria:

-1988-1991

-5000+ acres of burning

-500+ acres of high-severity

-Within selected ecoregions

->25% of selected forest types

# mtbs fire perimeters
mtbs_full <- st_read('data/mtbs_perims_DD.shp') %>% 
  st_transform(., crs=crs)
## Reading layer `mtbs_perims_DD' from data source 
##   `G:\Other computers\My Laptop\Documents\Grad School\Research\ConiferRegeneration\data\mtbs_perims_DD.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 29533 features and 22 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -166.1885 ymin: 17.94736 xmax: -65.33821 ymax: 70.15893
## Geodetic CRS:  NAD83
# filter fires by general region of interest and date range
mtbs_select <- mtbs_full %>% 
  mutate(state = str_sub(Event_ID,0,2),
         year = year(as.Date(Ig_Date))) %>% 
  filter(state %in% c("WA","ID","MT","WY","SD"),
         between(Ig_Date, as.Date('1988-01-1'), as.Date('1991-12-31'))) 

mapview(mtbs_select,zcol="year")
# function to group adjoining fire polygons to ensure patches are contiguous
group_fires <- function(mtbs_year) {

  # join the polygons with themselves, and remove those that do not join with any besides themselves
  combined<- st_join(mtbs_year, mtbs_year, join=st_is_within_distance, dist = 180, left = TRUE,remove_self = TRUE) %>% 
    drop_na(Event_ID.y)%>% 
    dplyr::select(Event_ID.x,Event_ID.y)
  
  if(nrow(combined)>=1){
    
    # partition data into that that has overlap, and that that does not
    overlap <- mtbs_year %>%
      filter(Event_ID %in% combined$Event_ID.x)
    no_overlap <- mtbs_year %>%
      filter(!(Event_ID %in% combined$Event_ID.x))
    
    print(paste0("there are ",nrow(overlap)," overlapping polygons"))
    
    # join all overlapping features, and buffer with radius to ensure proper grouping
    overlap_union <- st_union(overlap) %>%
    st_buffer(190)
    
    # break apart the joined polygons into their individual groups
    groups <- st_as_sf(st_cast(overlap_union ,to='POLYGON',group_or_split=TRUE)) %>%
      mutate(year = mean(mtbs_year$year),
             Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
      rename(geometry = x)
    
    print(paste0("polygons formed into ",nrow(groups)," groups"))
    
    # join back with original dataset to return to unbuffered geometry
    grouped_overlap <- st_join(overlap,groups,left=TRUE)
    
    # arrange by the new grouping
    joined_overlap_groups <- grouped_overlap %>%
      group_by(Fire_ID) %>%
      tally()%>%
      st_buffer(1) %>%
      dplyr::select(Fire_ID) %>%
      mutate(year = mean(mtbs_year$year))
    
    # add new ID to the freestanding polygons
    no_overlap_groups <- no_overlap %>%
      mutate(Fire_ID = str_c("Fire_",nrow(groups)+c(1:nrow(no_overlap)),"_",year)) %>%
      dplyr::select(Fire_ID,year)
    
    # join the new grouped overlap and the polygons without overlap
    fires_export <- rbind(joined_overlap_groups,no_overlap_groups)
    return(fires_export)
    
    } else { 
      
      print("no overlapping polygons")
      
      fires_export <- mtbs_year %>%
        mutate(Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
        dplyr::select(Fire_ID,year)
      
      return(fires_export)
  }
}
# group adjacent polygons within each fire year
fires_88 <- group_fires(mtbs_select %>%  filter(year == 1988))
## [1] "there are 22 overlapping polygons"
## [1] "polygons formed into 7 groups"
fires_89 <- group_fires(mtbs_select %>%  filter(year == 1989))
## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_90 <- group_fires(mtbs_select %>%  filter(year == 1990))
## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_91 <- group_fires(mtbs_select %>%  filter(year == 1991))
## [1] "no overlapping polygons"
# join each fire year, filter by area
mtbs_grouped <- rbind(fires_88,fires_89,fires_90,fires_91)%>%
  mutate(aream2 = as.numeric(st_area(geometry)),
         areaha = aream2/10000,
         areaac = areaha*2.471)%>% 
  filter(areaac >5000)
# assign proportions of forest type to each fire polygon
fires_join <- st_join(mtbs_grouped,eco_select,join=st_intersects,left=FALSE,largest=TRUE) %>% 
  left_join(., exact_extract(conus_forestgroup,mtbs_grouped, append_cols = TRUE, max_cells_in_memory = 3e+08, 
                             fun = function(value, coverage_fraction) {
                               data.frame(value = value,
                                          frac = coverage_fraction / sum(coverage_fraction)) %>%
                                 group_by(value) %>%
                                 summarize(freq = sum(frac), .groups = 'drop') %>%
                                 pivot_wider(names_from = 'value',
                                             names_prefix = 'freq_',
                                             values_from = 'freq')}) %>%
              mutate(across(starts_with('freq'), replace_na, 0)))
 
# remove unnecessary columns, cleanup names
# filter to remove majority other forest types, majority not-forested cover
fires <- fires_join %>% 
  dplyr::select("Fire_ID","year","areaha","areaac","US_L3NAME","freq_0","freq_200","freq_220","freq_260","freq_280") %>% 
  rename("fire_name"="Fire_ID",
         "fire_acres"="areaac",
         "ecoregion" = "US_L3NAME",
         "freq_df"="freq_200",
         "freq_pp"="freq_220",
         "freq_fs"="freq_260",
         "freq_lpp"="freq_280") %>% 
  mutate(freq_allother = 1-(freq_0 + freq_df+freq_pp+freq_fs+freq_lpp),
         freq_forested = 1- freq_0,
         freq_ideal = freq_df+freq_fs+freq_lpp)%>% 
  mutate(across(starts_with('freq'), round,2))%>% 
  filter(freq_ideal > 0.25)
# import all mtbs rasters via a list
rastlist <- list.files(path = "data/MTBSmosaic", pattern='.tif', all.files=TRUE, full.names=TRUE)
allrasters <- lapply(rastlist, raster)
names(allrasters) <- str_c("y", str_sub(rastlist,28,31))

# create empty dataframe
severity_list <- list()

# loop through mtbs mosasics for 1986-1991
# extract burn severity raster for all selected fires
# calculate burn severity percentages for each fire
for (i in names(allrasters)){
  mtbs_year <- allrasters[[i]]
  fire_year <- filter(fires, year==str_sub(i,2,5)) 
  raster_extract <- exact_extract(mtbs_year,fire_year, max_cells_in_memory = 3e+09,coverage_area=TRUE)
  names(raster_extract) <- fire_year$fire_name 
  
  output_select <- bind_rows(raster_extract, .id = "fire_name")%>%
    group_by(fire_name , value) %>%
    summarize(total_area = sum(coverage_area)) %>%
    group_by(fire_name) %>%
    mutate(proportion = total_area/sum(total_area))%>% 
    dplyr::select("fire_name","value","proportion") %>% 
    spread(.,key="value",value = "proportion")
  
  severity_list[[i]] <- output_select
}

# combine extracted raster datasets
severity_df <- do.call(rbind, severity_list) 

# join burn severity % to fires polygons
# fix naming
# filter dataset for 500 acres high severity
fires_severity <- left_join(fires,severity_df,by="fire_name")%>% 
  rename(noburn= "1",lowsev = "2", medsev = "3", highsev = "4",regrowth = "5", error = "6") %>% 
  dplyr::select(- "NaN",-"regrowth",-"error") %>% 
  mutate(highsev_acres = fire_acres*highsev)%>% 
  filter(highsev_acres > 500) 

fires_select <- fires_severity %>%
  left_join(.,exact_extract(conus_forestgroup,fires_severity, 'mode', append_cols = TRUE, max_cells_in_memory = 3e+08)) 

fires_select$mode <- as.factor(fires_select$mode)

fires_select <- fires_select %>% 
    mutate(fire_foresttype = case_when(mode==200 ~ "Douglas-Fir",
                                       mode==220 ~ "Ponderosa",
                                       mode==260 ~ "Fir-Spruce",
                                       mode==280 ~ "Lodegepole Pine",
                                       TRUE ~ "Other")) 

Fire Dataset

Selected fires by year

# plot
mapview(fires_select, zcol = "year")

Selected fires by majority forest type

# plot
mapview(fires_select, zcol = "fire_foresttype")

Export Data

Merge Fire Polys for Analysis

fires_export <- fires_select %>% 
  dplyr::select('fire_name','year','fire_acres','ecoregion','fire_foresttype','geometry') %>% 
  rename(Fire_ID = fire_name) %>% 
  mutate(year = as.integer(year)) %>% 
  st_transform(., crs="EPSG:4326")

str(fires_export)
## sf [37 × 6] (S3: sf/tbl_df/tbl/data.frame)
##  $ Fire_ID        : chr [1:37] "Fire_2_1988" "Fire_3_1988" "Fire_4_1988" "Fire_5_1988" ...
##  $ year           : int [1:37] 1988 1988 1988 1988 1988 1988 1988 1988 1988 1988 ...
##  $ fire_acres     : num [1:37] 342005 777690 448911 5651 11945 ...
##  $ ecoregion      : chr [1:37] "Middle Rockies" "Middle Rockies" "Middle Rockies" "Idaho Batholith" ...
##  $ fire_foresttype: chr [1:37] "Fir-Spruce" "Lodegepole Pine" "Lodegepole Pine" "Douglas-Fir" ...
##  $ geometry       :sfc_GEOMETRY of length 37; first list element: List of 4
##   ..$ :List of 49
##   .. ..$ : num [1:40676, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:1047, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:1164, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:6, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:171, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:880, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:998, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110 -110 -110 -110 44.7 ...
##   .. ..$ : num [1:55, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:1132, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:5, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110 -110 -110 -110 44.8 ...
##   .. ..$ : num [1:80, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110 -110 -110 -110 44.7 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.8 ...
##   .. ..$ : num [1:950, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:6, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:811, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:5, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110 -110 -110 -110 44.6 ...
##   .. ..$ : num [1:1558, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.8 ...
##   .. ..$ : num [1:111, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.7 ...
##   .. ..$ : num [1:598, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:34, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:267, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:135, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:19, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.7 ...
##   .. ..$ : num [1:62, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.7 ...
##   .. ..$ : num [1:53, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.8 ...
##   .. ..$ : num [1:5, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:44, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.8 ...
##   .. ..$ : num [1:68, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:26, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:20, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.8 ...
##   .. ..$ : num [1:5, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.7 ...
##   .. ..$ : num [1:204, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.7 ...
##   .. ..$ : num [1:44, 1:2] -110 -110 -110 -110 -110 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.7 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.7 ...
##   .. ..$ : num [1:4, 1:2] -110.1 -110.1 -110.1 -110.1 44.7 ...
##   ..$ :List of 1
##   .. ..$ : num [1:5308, 1:2] -110 -110 -110 -110 -110 ...
##   ..$ :List of 1
##   .. ..$ : num [1:140, 1:2] -110 -110 -110 -110 -110 ...
##   ..$ :List of 1
##   .. ..$ : num [1:169, 1:2] -110 -110 -110 -110 -110 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:5] "Fire_ID" "year" "fire_acres" "ecoregion" ...

Export

# st_write(fires_export, "output/", "fires_export.shp", driver = 'ESRI Shapefile')